# install.packages("devtools")
# devtools::install_github("MathiasHarrer/dmetar")

library(dmetar)
library("meta")
library("metafor")
library("PerformanceAnalytics")

GLPT1DMMETA_r= read.rm5(
  "GLP-1 for Type 1 Diabetes- metareg.rm5",
  sep = ",",
  quote = "\"",
  "GLP-1 for Type 1 Diabetes- metareg",
  numbers.in.labels = TRUE,
  debug = 0
)

c2o1 <- metacr(GLPT1DMMETA_r,comp.no=2, outcome.no=1)
year21 <- c(2016, 2016, 2016, 2019, 2018, 2016, 2019, 2020, 2018, 2015, 2020, 2012, 2016, 2016, 2016, 2016, 2016, 2016, 2021, 2022)
#Location:   MN = 1, NA = 2, Europe = 3, Asia = 4
loc21 <- c(1, 1, 1, 3, 2, 3, 3, 3, 2, 3, 1, 4, 2, 2, 2, 1, 1, 1, 1, 3)
#Design:  Parallel =1 , Crossover = 2
des21 <- c(1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2)
#Dose: 
dose21 <- c(0.6, 1.2, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.2, 1.8, 0.9, 0.6, 1.2, 1.8, 0.6, 1.2, 1.8, 1.8, 1.2)
#Risk Low = 1, Moderate/Uncertain = 2, High = 3
rob21 <- c(1, 1, 1, 2, 2, 1, 2, 1, 2, 2, 3, 3, 2, 2, 2, 1, 1, 1, 2, 2)
#cpep: only new onset or c-peptide positive: 1, others: 2
cpep21 <- c(2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1)
age21 <- c(43.2, 43.2, 43.2, 50.4, 46.7, 48, 27, 46.5, 35.8, 37.8, 46, 48.5, 44, 44, 44, 43.7, 43.7, 43.7, 29.5, 33.6) 
bmi21 <- c(28.9, 28.9, 28.9, 29, 28.9, 30.1, 24, 29.5, 30.5, 23.4, 31.6, 22.8, 29, 29, 29, 29.5, 29.5, 29.5, 24.1, 22.9)
a1c21 <- c(8.08, 8.08, 8.08, 8.2, 7.82, 8.7, 8.4, 8.15, 7.4, 8.75, 7.89, 7.65, 7.57, 7.57, 7.57, 8.16, 8.16, 8.16, 7.25, 6.8)

factors21 <- array(c(year21, loc21, des21, dose21, rob21, cpep21, age21, bmi21, a1c21), dim=c(20,9))
chart.Correlation(factors21)

c2o1d= data.frame(TE= c(c2o1$TE), seTE= c(c2o1$seTE), 
                  year= c(year21),
                  loc= c(loc21),
                  des= c(des21),
                  dose= c(dose21),
                  rob= c(rob21),
                  cpep= c(cpep21),
                  age= c(age21),
                  bmi= c(bmi21),
                  a1c= c(a1c21))

c2o1.mi <- multimodel.inference(TE= "TE", seTE = "seTE", data = c2o1d,
                     predictors = c("loc", "rob", "dose", "des",
                                    "cpep","bmi","a1c","year","age"),
                     interaction=FALSE)
c2o1.mi
#
c2o1.mreg <- rma(yi = TE,
                sei = seTE,
                data = c2o1,
                method = "ML",
                mods = ~year21+rob21,
                test = "knha")
c2o1.mreg
c2o1.mreg1 <- rma(yi = TE,
                 sei = seTE,
                 data = c2o1,
                 method = "ML",
                 mods = ~year21,
                 test = "knha")

c2o1.mreg1
anova(c2o1.mreg, c2o1.mreg1)



c2o2 <- metacr(GLPT1DMMETA_r,comp.no=2, outcome.no=2)
year22 <- c(2016, 2016, 2016, 2019, 2018, 2016, 2020, 2018, 2015, 2020, 2016, 2016, 2016, 2016, 2016, 2016, 2021, 2022)
#Location:   MN = 1, NA = 2, Europe = 3, Asia = 4
loc22 <- c(1, 1, 1, 3, 2, 3, 3, 2, 3, 1, 2, 2, 2, 1, 1, 1, 1, 3)
#Design:  Parallel =1 , Crossover = 2
des22 <- c(1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2)
#Dose: 
dose22 <- c(0.6, 1.2, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.2, 1.8,  0.6, 1.2, 1.8, 0.6, 1.2, 1.8, 1.8, 1.2)
#Risk Low = 1, Moderate/Uncertain = 2, High = 3
rob22 <- c(1, 1, 1, 2, 2, 1, 1, 2, 2, 3, 2, 2, 2, 1, 1, 1, 2, 2)
#cpep: only new onset or c-peptide positive: 1, others: 2
cpep22 <- c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1)
age22 <- c(43.2, 43.2, 43.2, 50.4, 46.7, 48, 46.5, 35.8, 37.8, 46, 44, 44, 44, 43.7, 43.7, 43.7, 29.5, 33.6) 
bmi22 <- c(28.9, 28.9, 28.9, 29, 28.9, 30.1, 29.5, 30.5, 23.4, 31.6, 29, 29, 29, 29.5, 29.5, 29.5, 24.1, 22.9)
a1c22 <- c(8.08, 8.08, 8.08, 8.2, 7.82, 8.7, 8.15, 7.4, 8.75, 7.89, 7.57, 7.57, 7.57, 8.16, 8.16, 8.16, 7.25, 6.8)
# 
# factors22 <- array(c(year22, loc22, des22, dose22, rob22, cpep22), dim=c(18,5))
# chart.Correlation(factors22)

c2o2d= data.frame(TE= c(c2o2$TE), seTE= c(c2o2$seTE), 
                  year= c(year22),
                  loc= c(loc22),
                  des= c(des22),
                  dose= c(dose22),
                  rob= c(rob22),
                  cpep= c(cpep22),
                  age= c(age22),
                  bmi= c(bmi22),
                  a1c= c(a1c22))

c2o2.mi <- multimodel.inference(TE= "TE", seTE = "seTE", data = c2o2d,
                                predictors = c("loc", "rob", "dose", "des",
                                               "cpep","bmi","a1c","year","age"),
                                interaction=TRUE)
c2o2.mi
#
c2o2.mreg <- rma(yi = TE,
                 sei = seTE,
                 data = c2o2,
                 method = "ML",
                 mods = ~dose22+loc22+rob22,
                 test = "knha")
c2o2.mreg
c2o2.mreg1 <- rma(yi = TE,
                  sei = seTE,
                  data = c2o2,
                  method = "ML",
                  mods = ~dose22+loc22+rob22+cpep22,
                  test = "knha")

c2o2.mreg1
anova(c2o2.mreg, c2o2.mreg1)




c2o3 <- metacr(GLPT1DMMETA_r,comp.no=2, outcome.no=3)
year23 <- c(2016, 2016, 2016, 2019, 2016, 2019, 2020, 2018, 2015, 2020, 2012, 2016, 2016, 2016, 2016, 2016, 2016, 2022)
#Location:   MN = 1, NA = 2, Europe = 3, Asia = 4
loc23 <- c(1, 1, 1, 3, 3, 3, 3, 2, 3, 1, 4, 2, 2, 2, 1, 1, 1, 3)
#Design:  Parallel =1 , Crossover = 2
des23 <- c(1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2)
#Dose: 
dose23 <- c(0.6, 1.2, 1.8, 1.8,  1.8, 1.8, 1.8, 1.8, 1.2, 1.8, 0.9, 0.6, 1.2, 1.8, 0.6, 1.2, 1.8, 1.2)
#Risk Low = 1, Moderate/Uncertain = 2, High = 3
rob23 <- c(1, 1, 1, 2, 1, 2, 1, 2, 2, 3, 3, 2, 2, 2, 1, 1, 1, 2)
#cpep: only new onset or c-peptide positive: 1, others: 2
cpep23 <- c(2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1)
age23 <- c(43.2, 43.2, 43.2, 50.4, 48, 27, 46.5, 35.8, 37.8, 46, 48.5, 44, 44, 44, 43.7, 43.7, 43.7, 33.6) 
bmi23 <- c(28.9, 28.9, 28.9, 29, 30.1, 24, 29.5, 30.5, 23.4, 31.6, 22.8, 29, 29, 29, 29.5, 29.5, 29.5, 22.9)
a1c23 <- c(8.08, 8.08, 8.08, 8.2, 8.7, 8.4, 8.15, 7.4, 8.75, 7.89, 7.65, 7.57, 7.57, 7.57, 8.16, 8.16, 8.16, 6.8)

# factors23 <- array(c(year23, loc23, des23, dose23, rob23, a1c23), dim=c(18,6))
# chart.Correlation(factors23)

c2o3d= data.frame(TE= c(c2o3$TE), seTE= c(c2o3$seTE), 
                  year= c(year23),
                  loc= c(loc23),
                  des= c(des23),
                  dose= c(dose23),
                  rob= c(rob23),
                  cpep= c(cpep23),
                  age= c(age23),
                  bmi= c(bmi23),
                  a1c= c(a1c23))

c2o3.mi <- multimodel.inference(TE= "TE", seTE = "seTE", data = c2o3d,
                                predictors = c("loc", "rob", "dose", "des",
                                               "cpep","bmi","a1c","year","age"),
                                interaction=FALSE)
c2o3.mi
#
c2o3.mreg <- rma(yi = TE,
                 sei = seTE,
                 data = c2o3,
                 method = "ML",
                 mods = ~cpep23+dose23+des23+loc23+a1c23,
                 test = "knha")
c2o3.mreg
c2o3.mreg1 <- rma(yi = TE,
                  sei = seTE,
                  data = c2o3,
                  method = "ML",
                  mods = ~dose23,
                  test = "knha")

c2o3.mreg1
anova(c2o3.mreg, c2o3.mreg1)





c2o4 <- metacr(GLPT1DMMETA_r,comp.no=2, outcome.no=4)
year24 <- c(2016, 2016, 2016, 2019, 2016, 2020, 2018, 2015, 2020, 2016, 2016, 2016, 2016, 2016, 2016, 2021, 2022)
#Location:   MN = 1, NA = 2, Europe = 3, Asia = 4
loc24 <- c(1, 1, 1, 3, 3, 3, 2, 3, 1, 2, 2, 2, 1, 1, 1, 1, 3)
#Design:  Parallel =1 , Crossover = 2
des24 <- c(1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2)
#Dose: 
dose24 <- c(0.6, 1.2, 1.8, 1.8, 1.8, 1.8, 1.8, 1.2, 1.8,  0.6, 1.2, 1.8, 0.6, 1.2, 1.8, 1.8, 1.2)
#Risk Low = 1, Moderate/Uncertain = 2, High = 3
rob24 <- c(1, 1, 1, 2, 1, 1, 2, 2, 3, 2, 2, 2, 1, 1, 1, 2, 2)
#cpep: only new onset or c-peptide positive: 1, others: 2
cpep24 <- c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1)
age24 <- c(43.2, 43.2, 43.2, 50.4, 48, 46.5, 35.8, 37.8, 46, 44, 44, 44, 43.7, 43.7, 43.7, 29.5, 33.6) 
bmi24 <- c(28.9, 28.9, 28.9, 29, 30.1, 29.5, 30.5, 23.4, 31.6, 29, 29, 29, 29.5, 29.5, 29.5, 24.1, 22.9)
a1c24 <- c(8.08, 8.08, 8.08, 8.2, 8.7, 8.15, 7.4, 8.75, 7.89, 7.57, 7.57, 7.57, 8.16, 8.16, 8.16, 7.25, 6.8)
# 
# factors24 <- array(c(year24, loc24, des24, dose24, rob24, a1c24), dim=c(17,6))
# chart.Correlation(factors24)

c2o4d= data.frame(TE= c(c2o4$TE), seTE= c(c2o4$seTE), 
                  year= c(year24),
                  loc= c(loc24),
                  des= c(des24),
                  dose= c(dose24),
                  rob= c(rob24),
                  cpep= c(cpep24),
                  age= c(age24),
                  bmi= c(bmi24),
                  a1c= c(a1c24))

c2o4.mi <- multimodel.inference(TE= "TE", seTE = "seTE", data = c2o4d,
                                predictors = c("loc", "rob", "dose", "des",
                                               "cpep","bmi","a1c","year","age"),
                                interaction=FALSE)
c2o4.mi
#
c2o4.mreg <- rma(yi = TE,
                 sei = seTE,
                 data = c2o4,
                 method = "ML",
                 mods = ~dose24+a1c24,
                 test = "knha")
c2o4.mreg
c2o4.mreg1 <- rma(yi = TE,
                  sei = seTE,
                  data = c2o4,
                  method = "ML",
                  mods = ~dose24,
                  test = "knha")

c2o4.mreg1
anova(c2o4.mreg, c2o4.mreg1)
